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Abstract — In considering tlie reliability of numerical pro- 
grams, it is normal to "limit our study to the semantics 
dealing with numerical precision" (Martel, 2005). On the other 
hand, there is a great deal of work on the reliability of 
programs that essentially ignores the numerics. The thesis 
of this paper is that there is a class of problems that fall 
between these two, which could be described as "does the low- 
level arithmetic implement the high-level mathematics". Many 
of these problems arise because mathematics, particularly 
the mathematics of the complex numbers, is more difficult 
than expected: for example the complex function log is not 
continuous, writing down a program to compute an inverse 
function is more complicated than just solving an equation, 
and many algebraic simplification rules are not universally 
valid. 

The good news is that these problems are theoretically 
capable of being solved, and are practically close to being solved, 
but not yet solved, in several real-world examples. However, 
there is still a long way to go before implementations match 
the theoretical possibilities. 

Aeyworefa -Verification; algebra; branch cut; singularity; 

I. Introduction 

It is customary, even though not often explicitly stated, to 
think that programming errors in numerical programs come 
in three distinct flavours, which we can categorise as follows. 

1. Blunder (of the coding variety) This is the sort 
of error traditionally addressed in "program verifi- 
cation": are all array elementsproperly initialised 
before use, "fence post" errorsl^ are array bounds 
always respected etc.? These problems are essen- 
tially independent of the numerics of the problem, 
and indeed are normally taught/described in integer 
contexts. 

2. Parallelism Issues of deadlocks or races occur- 
ring due to the parallelisation of an otherwise 
correct sequential program. Again, these problems 
are essentially independent of the numerics of the 
problem. 

'From the old puzzle "A farmer wishes to make a lOO-metre fence with 
supporting posts every 10 metres — how many posts are needed", to which 
the answer is 1 1 , since each end needs a post. 



3. Numerical Do truncation and round-off errors, 
individually or combined, mean that the program 
computes approximations to the "true" answers 
which are out of tolerance. This is the area tra- 
ditionally addressed in Numerical Analysis. There 
are really two subquestions here: the rounding 
question, i.e. does Hjeee (or whatever other arith- 
metic we are using) approximate R sufficiently 
well, and the truncation error question, e.g. is 
the discretisation h small enough that it is the 
mathematical e or is equivalent to X^-T- Un- 
fortunately the two interact; for example reducing 
h in f'{x) w fi^+^y-t{^') [Q reduce the truncation 
error increases the rounding problems. 

We note that [l], taken as a specimen of the verification 
literature, contains 30 papers, of which only [2 1 deals with 
strictly numerical issues, four with parallelism issues, and 
the rest (83%) with the first kind. 

It is the thesis of this paper that there is a fourth kind of 
error, which we can describe as follows 

4. Manipulation A piece of algebra, which is "obvi- 
ously correct", turns out not to be correct when 
interpreted, not as abstract algebra, but as the 
manipulation of functions R — ^ R or C — > C. 

Note: throughout this paper we take the standard definitions 
of the branch cuts of the elementary functions from |3 |, |4J 
(as tightened in Jsj). Other definitions would have different, 
but not fewer, problems. We also use the Anglo-Saxon con- 
vention that log etc. (and ^) denote single-valued functions 
(log 1 = 0, \/4 = 2), whereas Log etc. denote multi-valued 
functions (Log(l) = {2kTTi : fc G Z}, Sqrt(4) = {2, -2}). 

The problems we are going to describe arise largely 
(though not entireljQ) from complex numbers, and it is 
sometimes said "real programs don't use complex numbers", 
despite the fact that the introduction of COMPLEX into 
Fortran II was probably the first instance of a language data 
type that did not correspond to a machine data type. The 
authors know of several major uses of complex numbers 
and complex analysis, in particular many problems which 

^See section ITvl for a counter-example. 



arise in fluid dynamics, where two-dimensional real space 
= {{x,y)} is viewed as the complex plane C — {z = 
X + iy}. It is then normal to map this copy of C to another 
(in which the variable is traditionally denoted w or Q where 
the problem is easier to solve. Such an analytic map z ^ w 
is termed a conformal map. 

II. The Kahan Example 

This example comes from ||6] pp. 187-189], and the 
ultimate motivation is fluid flow in a slotted strip (2 space), 
which we wish to transform to a more convenient shape (w 
space). 

With the usual definitions, the necessary conformal map 
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w = g[z) := 2 arccosn 1 + -5- I ^ arccosn 



is only the same as the ostensibly more efficient 

w^q{z) := 2arccosh ( 2(z + 8)4 
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if we avoid the teardrop shaped area 

2 = x+iy : \y\ < 
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as shown by Figure [U 
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developed at Bath and to be included with Maple 17 does 
isolate the curve 



y^± 



(a; + 3)2(-2a;- 9) 



2a; + 5 



(4) 



with the appropriate a; range as a potential obstacle (it is 
the branch cut of q). However, it is just one of a set of cuts 
returned by the code. The plotting option in the package 
produces Figure |2] which presents the teardrop and the entire 
real axis as potential cuts. 




Figure 2. Plot of the potential branch cuts for g{z) = q{z) produced by 
the BranchCuts package. 




Figure 1. Plots of the real and imaginary parts of g{z) — q{z). 

In fact, most computer algebra systems will refuse, these 
days, to convert one into the other, but this does not 
constitute a proof of difference. 

Challenge 1: Demonstrate automatically that g and q are 
not equal, by producing a z at which they give different 
results. 

The technology described in fT] would answer this question 
by decomposing the complex plane into regions, via means 
of cylindrical algebraic decomposition (CAD) with respect 
to the branch cuts of the functions, and testing the identity 
at a sample point in each region 

A fully-automated implementation for this example has 
yet to be produced since the geometry can be quite involved. 
(See section HI] for details.) However, progress is currently 
being made on the problem. The BranchCuts package fS] 



The package calculates cuts for functions by mapping the 
defining branch cuts of a function by the argument. The 
output is generated using Maple's solve facilities, and 
the user can choose to look for solutions in the complex 
variable, or to first split into real and imaginary parts and 
work over the reals. The first method is quicker and more 
widely applicable, but the second produces more intuitive 
output including the algebraic description of the teardrop in 
equation (HJl. 

The package identifies the potential branch cuts of a 
composition of functions (a sum, a product or an equality) 
as the union of the cuts for individual components. Hence 
the identification of the real axis as a potential obstacle is 
not surprising since the individual terms do have branch cuts 
here, with the resulting discontinuities happening to cancel 
out in the composition. 

However, the output described here would not be suitable 
for the technology of [7 1 since the input into CAD must be 
semi-algebraic. We can modify the code to just return the 
polynomial equalities and inequalities that define each set of 
cuts. For this example, there are 7 such sets, one of which 
consists of the three relations below. 

4:y{2y^x + 2x^ + 5y^ + 21x^ + 72a; + 81) = 



+9y2 - 63a;2 - 135a; - 108) < 
-63y^ + 22hx^ + 324x < 



(5) 



Figure |3] gives a plot of these three curves. By testing sample 
points we see where the inequalities are satisfied and infer 
that the branch cut defined is the teardrop along with the 
real axis above —3. These issues and the implementation of 
the BranchCuts package are discussed further in |,8J. 




Figure 3. Plot of the information in (|5]. The solid line is the equality, the 
dashed line the first inequality and the dotted line the second. 

If we pass the full set of polynomials to CAD (ignoring 
whether they are equalities or inequalities) then clearly 
a lot more information will be processed than required. 
Using the variable ordering y > x and the command 
CylindricalAlgebraicDecompose within Maple 16 
||9J, this produces a CAD of 409 cells. Given (O we might 
hope for a minimal CAD of 13 cells, or if we accept that 
the real axis must be included in any calculations then a 
minimal CAD of 19 cells, (see Figure |4|i. 
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Figure 4. Possible minimal CADs for Kahan's example. The dots indicate 
sample points for a cell. 

Note that there is a simplification of g valid over the whole 
complex plane : g can legitimately be rewritten to 

, oi /l V37TT2(VIT3 + yif \ 



The technology in Q can show this, i.e. Vz G Cg{z) = 
h{z), but there are multiple square roots requiring denesting 
1261 §4.3] and (as formulated) square roots in the denomina- 
tor. Indeed the standard tools of Maple 16 currently get this 
wrong: coulditbe (goh) ; returns true, which ought 
to indicate that there is a counter-example. 

Challenge 2: Demonstrate automatically that g and h are 
equal. 

Although the technology in l?], implemented in a mixture 
of Maple and QEPCAD (though we are working on a purely 
Maple implementation based on |9|), may be able to meet 
this challenge in time, we would be left with the problem 
of trusting the underlying demonstration code. So there is 
the additional problem of translating this methodology into 
a tool such as MetiTarski [lO]. 

III. JouKowsKi Examples 
Consider the Joukowski map fTTi pp. 294-298]: 



/ 



1 

'^2 



1 

z+- 
z 



(7) 



A. Joukowski (a) 

Lemma 1: f is injective as a function from D := {z : 
\z\ > 1}. 

If 2; C then 1/z i-^ and there are no other pre-images 
of (since the algebraic inverse of (|7) is the solution of a 
quadratic). If \z\ > 1, then \l/z\ < 1, so z is unique in D. 

In fact / is a bijection from Z) to C* := C \ [—1, 1], and 
hence has an inverse. 

Of course, (|7]i is the conformal map C C that equates 
to the map 



/fl : {x, y) ^ 
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^2 2 x2 +y2' 2 " 2 x^ +y^ 
R2 — > R2. However, it is not obvious from (|8]l alone that 



Jr is a bijection D 

Va;iVa;2VyiVy2 



i.e. that 



xl + yl > 1 A xj + yj > lA 
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Challenge 3: Demonstrate automatically the truth of (|9]l, 
which is also L12, (49)]. 

We have been unable to do this with either the QEPCAD 
fT3\ implementation of Partial Cylindrical Algebraic De- 
composition 1 14 1 or the Maple implementation of Cylindri- 
cal Algebraic Decomposition via triangular decomposition 
0. 

Brown ifTSl suggested writing problems of the form 
Vv.4 (/ = A g = 0) as -<3vA A = A .g = 0) . 
This splits into two clauses: fT2] §7.2], which can be seen 
to have equational constraints lil^ . 



Figure 5. Maple's solve on inverting Joukowski 
> [solve(zeta = l/2*(z+l/z), z) ] ; 



1 



Figure 6. Maple's actual solve on inverting injective Joukowski 

> [solve(zeta = l/2*(z+l/z), z)]\ 
assuming abs(z) > 1 



1 



Even these are beyond QEPCAD and Maple currently. 
However, Brown |15 | has been able to recast these clauses 
(manually) to make them amenable to QEPCAD, and indeed 
solved the problem in under 12 seconds. 

Challenge 4: Automate these techniques and transforms. 
Having established (or not) that / is a bijection D — > C^, 
we want its inverse. Formally, this is trivial, as one referee 
said 

The inverse of Joukowski is the solution of a 
quadratic with the usual sign ambiguity: 



if C = 5 (z + 7), then 2zC = + 1 and z = C± v^C^ - 1. 
This is easily within the grasp of computer algebra, as seen 
in Figure |5] The only challenge might be the choice implicit 
in the ± symbol: which do we choose? 

Unfortunately, the answer is "neither", or at least "neither, 
uniformly". The true answer is that, for / a bijection from 
{z : \z\ > 1} to C \ [—1,1], its inverse is 



/i(C) = C 



+v e^ 5(c) >o 

+ VCZl 3(C) = A 5R(C) > 1 
-VC^ 5(C) = 0A5R(C)<-1 



(10) 



In fact, a better (at least, free of case distinctions) definition 
is 

' ' (11) 



/2(C) = C + v/C^v/cTT. 



The techniques of Q are able to verify ( fTTT l. in the sense 
of showing that /2(/(z)) — z is the zero function on {z : 
|.|>1}. 

Challenge 5: Derive automatically, and demonstrate the 
validity of, either (fTOl i or (flU . In terms of Maple, we would 
want to see Figure [T] rather than the actual Figure |6] 



Figure 7. Ideal Maple solve on inverting injective Joukowski 

> solve(zeta = l/2*(z+l/z), z)\ 
assuming abs(z) > 1 



In terms of derivation, the techniques of |17 | seem worthy 
of investigation, but the first author has been unable to do 
this derivation satisfactorily by this route. 

B. Joukowski (b) 

Here the function is again given by O. 
Lemma 2: f is injective as a function from H := {z : 
5z > 0}. 

As in Lemma [T] if z H> C then 1/z H> C> ™d there are no 
other pre-images of C- If S(z) > 0, 3(l/z) < 0, and / in 
therefore injective from H. 

In fact, / is a bijection from H to C\{{—oo, — 1]U[1, 00)), 
and hence has an inverse. 

Again, it is not obvious from (|8]l alone that fjj is a 
bijection, now from {{x,y)\y > 0}, i.e. that 
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(xi = X2 A J/1 = y2)- 

Challenge 6: Demonstrate automatically the truth of (fT2] l. 
It is likely that the ideas of ifTSl can do this, but again these 
need automation. 

We have the same challenge over the inverse of /: again 
formally it is f^^—(^± y^(^~~T, and the only challenge is 
the ± symbol: which do we choose? Here ifTTI (5.1-13), p. 
298] argues for 



/3(C) = C+ VC^ VC+l- (13) 

argG(-7r/2,7r/2] argg(0,7r] 

Challenge 7: Find a way to represent functions such as 

TcTi 

arge(0.7r] 

Fortunately this one is soluble in this casqj, we can write 
— i -\/— C — 1 , and the latter is the normal 

arge(0,7r] argG(-Tr/2,Tr/2] 

sqrt function of [3'1. Hence we have an inverse function 



/4(C) = c + vc^^V-c-i. 



(14) 



Challenge 8: Demonstrate automatically that this is an 
inverse to / on {z : 5z > 0}. 

There is also the problem of deducing (fl4] l. One could try 
automatic deduction on the lines of [ 17 1, but there is another 
possibility: heuristic search followed by verification the 
verification were feasible ifTSI . 



C + VC^VC+T 



'And is probably soluble more generally, but the authors know of no 
general work on "alternative formulations". 



IV. A Real Example 

Just in case the reader thinks that the real numbers are 
immune from these problems, consider the addition rule for 
the inverse tangent, quoted as ( Q (4.4.34)] lH (4.24.15)]) 



Arctan(x) ± Arctan(?/) = Arctan 



x±y 



Despite the caveat in |4| that "The above equations are 
interpreted in the sense that every value of the left-hand 
side is a value of the right-hand side and vice versa", it is 
in fact the case that the 'obvious' two equations are true 
separately, viz . 



Arctan(a;) + Arctan(j/) 
Arctan(a;) — Arctan(y) 
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(15) 



(16) 



Consider (fTsT i: This is valid for the multi-valued Arctan, 
but for the single-valued arctan only when |1 — xy\ < 1, 
due to a "branch cut at infinity" of arctan. Nevertheless, the 
single-valued version of ( fTSl l is often cited as true: see for 
example [19, (5.2.5)]. 

Over the reals, this is a non-challenge, the techniques of 
fi] do solve it easily, and produce a counterexample. 

V. So WHY ARE THESE CHALLENGES? 

A. Complex functions and branch cuts 

These are difficult subjects, which have tended to be 
brushed under the carpet. The first truly algorithmic ap- 
proach is ten years old ( [|20l . refined in JT]), and has various 
difficulties. 

1) At its core is the use of Cylindrical Algebraic De- 
composition of to find the connected components 
of C^/^ \ {branch cuts}. The complexity of this is 
doubly exponential in N: upper bound of ) ||2il 
and lower bounds of 2^'" (|22l, EH. While better 
algorithms for the connected components problem are 
in principle known ( [!241 is d'-'^^^''), we do not 
know of any accessible implementations. 
Furthermore, we are clearly limited to small values of 
N, at which point looking at 0{. . .) complexity is of 
limited use. We note that the cross-over point between 
2(w-i)/3 is at TV = 21. A more detailed 
comparison is given in [21] . Hence there is a need 
for practical research on low-iV CyUndrical Algebraic 
Decomposition. 

2) While the fundamental branch cut of log is simple 
enough, being {z = x + iy\y = A x < 0}, actual 
branch cuts are messier Part of the branch cut of (2) 
is 



2x'^ + 21x^ 



- 72x + 2x2/2 + 5?/2 + 81 = A 
other conditions, 



whose solution accounts for the curious expression in 
(O. While there has been some progress in manipulat- 
ing such images of half-lines (described in [25J, [.26] ). 
there is almost certainly more to be done. 

B. Injectivity 

Lemmas [T] and |2] might seem to be statements about 
complex functions of one variable, so why do we need to 
handle (or fail to handle) statements about four real variables 
to prove them? There are three, rather distinct, reasons for 
this. 

1) The statements require the | • | function (Lemma [T]) 
or the S function (Lemma |2]i, neither of which are 
complex analytic functions. Hence some recourse to 
real analysis (and therefore twice as many variables) 
seems inevitable, though it would be nice to have a 
more formal statement and proof of this. 

2) Equations (|9]l and (fT2l l are the direct translations of 
the basic definition of injectivity. In practice, certainly 
if we were looking at functions R R, we would 
want to use the fact that the function concerned was 
continuous. 

Challenge 9: Find a better formulation of injectivity 
questions R^ R^, making use of the properties of 
the functions concerned (certainly continuity, possibly 
rationality). 

3) While equations (|9]l and ( fT2] l are statements from the 
existential theory of the reals, and so the theoreti- 
cally more efficient algorithms quoted in [21] are in 
principle applicable, the more modem developments 
described in [27] do not seem to be directly applicable. 
However, we can transform then into a disjunction of 
statements to each of which the Weak Positivstellen- 
satz [27 , Theorem 1] is applicable. 

Challenge 10: Solve these problems using the tech- 
niques of ll27l . 

VI. Conclusions 
The aim of this paper has been to demonstrate that 
translating mathematical problems into programs may re- 
quire some algebraic manipulations whose accuracy is not 
as obvious as one might think, and whose verification is 
currently not as straightforward as we would like, despite 
the fact that their correctness is, in principle, decidable. A 
summary is given in Table U 

These are, largely, concrete challenges that, we hope, will 
spur practical advances in this domain. 
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Table I 

Current state of these challenges 
State 

Mathematically solved by (7), 
geometiy taxes cuiTent solvers. 
Mathematically solved by fj\, 
branch cut determination 
not yet implemented. 
Mathematically solved by |28 etc.], 
geometry defeats current solvers. 
Under development. 

Mathematically solved (7|, geometry defeats current 
solvers, and is probably significantly harder than 
in the previous problems. 

Unknown: probably straightforward research project. 

Unknown: research project. 

Unknown: project for the authors of 1271 . 
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